1 Model

Loss curves are a standard actuarial technique for helping insurance companies assess the amount of reserve capital they need to keep on hand to cover claims from a line of business. Claims made and reported for a given accounting period are tracked seperately over time so that more recent periods with less claim development can use the behaviour of policies from earlier periods of time to help predict what total claims will be.

Total claim amounts from a simple accounting period are usually laid out in a single line with each column showing the total claim amount after that period of time. Subsequent accounting periods have less development and so the data takes a triangular shape, hence the term ‘loss triangles’. Using previous patterns, data in the upper part of the triangle can be used to predict values in the unknown lower triangle.

The ChainLadder package provides functionality to generate and use these loss triangles.

In this case study, we take a related but different approach: we model the growth of the losses in each accounting period as an increasing function of time, and use the model to estimate the parameters which determine the shape and form of this growth. We also use the sampler to estimate the values of the “ultimate loss ratio”, i.e. the ratio of the total claims on an accounting period to the total premium received to write those policies.

1.1 Overview

We will work with two different functional forms for the growth behaviour of the loss curves: a ‘Weibull’ model and a ‘loglogistic’ model:

\[ g(t \, ; \; \theta, \omega) = \frac{t^\omega}{t^\omega + \theta^\omega} \;\; (\text{Weibull}) \]

\[ g(t \, ; \; \theta, \omega) = 1 - \exp\left(-\left(\frac{t}{\theta}\right)^\omega\right) \;\; (\text{Log-logistic}) \]

2 Load Data

We will first look at the car insurance loss data from casact.org

### File was downloaded from http://www.casact.org/research/reserve_data/ppauto_pos.csv

ppauto_tbl <- read_csv("data/ppauto_pos.csv", progress = FALSE)
## Parsed with column specification:
## cols(
##   GRCODE = col_integer(),
##   GRNAME = col_character(),
##   AccidentYear = col_integer(),
##   DevelopmentYear = col_integer(),
##   DevelopmentLag = col_integer(),
##   IncurLoss_B = col_integer(),
##   CumPaidLoss_B = col_integer(),
##   BulkLoss_B = col_integer(),
##   EarnedPremDIR_B = col_integer(),
##   EarnedPremCeded_B = col_integer(),
##   EarnedPremNet_B = col_integer(),
##   Single = col_integer(),
##   PostedReserve97_B = col_integer()
## )
names(ppauto_tbl) <- ppauto_tbl %>% names %>% tolower

head(ppauto_tbl)
## # A tibble: 6 × 13
##   grcode                  grname accidentyear developmentyear developmentlag
##    <int>                   <chr>        <int>           <int>          <int>
## 1     43 IDS Property Cas Ins Co         1988            1988              1
## 2     43 IDS Property Cas Ins Co         1988            1989              2
## 3     43 IDS Property Cas Ins Co         1988            1990              3
## 4     43 IDS Property Cas Ins Co         1988            1991              4
## 5     43 IDS Property Cas Ins Co         1988            1992              5
## 6     43 IDS Property Cas Ins Co         1988            1993              6
## # ... with 8 more variables: incurloss_b <int>, cumpaidloss_b <int>,
## #   bulkloss_b <int>, earnedpremdir_b <int>, earnedpremceded_b <int>,
## #   earnedpremnet_b <int>, single <int>, postedreserve97_b <int>

In terms of modeling the work, we want to ensure the data we work with is a snapshot in time. To ensure this, we filter out all data available after 1997 and use the remaining data as our dataset.

Once we have created fits and made predictions, we use the later data as a way to validate the model.

use_grcode <- c(43,353,388,620)

carrier_full_tbl <- ppauto_tbl %>%
    mutate(acc_year   = as.character(accidentyear)
          ,dev_year   = developmentyear
          ,dev_lag    = developmentlag
          ,premium    = earnedpremdir_b
          ,cum_loss   = cumpaidloss_b
          ,loss_ratio = cum_loss / premium) %>%
    select(grcode, acc_year, dev_year, dev_lag, premium, cum_loss, loss_ratio)

carrier_snapshot_tbl <- carrier_full_tbl %>%
    filter(grcode %in% use_grcode
          ,dev_year < 1998)

We are looking at four insurers with the GRCODEs above. Before we proceed with any analysis, we first plot the data, grouping the loss curves by accounting year and faceting by carrier.

ggplot(carrier_snapshot_tbl) +
    geom_line(aes(x = dev_lag, y = loss_ratio, colour = as.character(acc_year))
             ,size = 0.3) +
    expand_limits(y = c(0,1)) +
    facet_wrap(~grcode) +
    xlab('Development Time') +
    ylab('Loss Ratio') +
    ggtitle('Snapshot of Loss Curves for 10 Years of Loss Development'
           ,subtitle = 'Private Passenger Auto Insurance for Single Organisation') +
    guides(colour = guide_legend(title = 'Cohort Year'))

We look at the chain ladder of the data, rather than looking at the loss ratios we just look at the dollar amounts of the losses.

snapshot_tbl <- carrier_snapshot_tbl %>%
    filter(grcode %in% use_grcode[1])

snapshot_tbl %>%
    select(acc_year, dev_lag, premium, cum_loss) %>%
    spread(dev_lag, cum_loss) %>%
    print
## # A tibble: 10 × 12
##    acc_year premium   `1`   `2`   `3`   `4`   `5`   `6`   `7`   `8`   `9`  `10`
## *     <chr>   <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1      1988     957   133   333   431   570   615   615   615   614   614   614
## 2      1989    3695   934  1746  2365  2579  2763  2966  2940  2978  2978    NA
## 3      1990    6138  2030  4864  6880  8087  8595  8743  8763  8762    NA    NA
## 4      1991   17533  4537 11527 15123 16656 17321 18076 18308    NA    NA    NA
## 5      1992   29341  7564 16061 22465 25204 26517 27124    NA    NA    NA    NA
## 6      1993   37194  8343 19900 26732 30079 31249    NA    NA    NA    NA    NA
## 7      1994   46095 12565 26922 33867 38338    NA    NA    NA    NA    NA    NA
## 8      1995   51512 13437 26012 31677    NA    NA    NA    NA    NA    NA    NA
## 9      1996   52481 12604 23446    NA    NA    NA    NA    NA    NA    NA    NA
## 10     1997   56978 12292    NA    NA    NA    NA    NA    NA    NA    NA    NA

We also look at the loss ratios in a similar fashion

snapshot_tbl %>%
    select(acc_year, dev_lag, premium, loss_ratio) %>%
    spread(dev_lag, loss_ratio) %>%
    print(digits = 2)
## # A tibble: 10 × 12
##    acc_year premium      `1`      `2`      `3`      `4`      `5`      `6`
## *     <chr>   <int>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1      1988     957 0.138976 0.347962 0.450366 0.595611 0.642633 0.642633
## 2      1989    3695 0.252774 0.472530 0.640054 0.697970 0.747767 0.802706
## 3      1990    6138 0.330727 0.792441 1.120886 1.317530 1.400293 1.424405
## 4      1991   17533 0.258769 0.657446 0.862545 0.949980 0.987909 1.030970
## 5      1992   29341 0.257796 0.547391 0.765652 0.859003 0.903752 0.924440
## 6      1993   37194 0.224310 0.535033 0.718718 0.808706 0.840162       NA
## 7      1994   46095 0.272589 0.584055 0.734722 0.831717       NA       NA
## 8      1995   51512 0.260852 0.504970 0.614944       NA       NA       NA
## 9      1996   52481 0.240163 0.446752       NA       NA       NA       NA
## 10     1997   56978 0.215732       NA       NA       NA       NA       NA
## # ... with 4 more variables: `7` <dbl>, `8` <dbl>, `9` <dbl>, `10` <dbl>

2.1 Quick Data Exploration

We are working with the loss ratio, so we recreate the chain ladder format but look at loss ratios instead of dollar losses.

ggplot(snapshot_tbl) +
    geom_line(aes(x = dev_lag, y = loss_ratio, colour = acc_year)
             ,size = 0.3) +
    expand_limits(y = 0) +
    xlab('Development Time') +
    ylab('Loss Ratio') +
    guides(colour = guide_legend(title = 'Cohort Year'))

3 Stan Model

3.1 Configure Data

We want to only use the data at a given snapshot, so we choose all data current to 1998. Thus, we have 10 years of development for our first 1988 cohort, and one less for each subsequent year. Our final cohort for 1997 has only a single year of development

cohort_tbl <- snapshot_tbl %>%
    select(acc_year, premium) %>%
    unique() %>%
    mutate(cohort_id = 1:n())

usedata_tbl <- snapshot_tbl

cohort_maxtime <- usedata_tbl %>%
    group_by(acc_year) %>%
    summarise(maxtime = max(dev_lag)) %>%
    arrange(acc_year) %>%
    .[[2]]

cohort_premium <- usedata_tbl %>%
    group_by(acc_year) %>%
    summarise(premium = unique(premium)) %>%
    .[[2]]

t_values <- usedata_tbl %>%
    select(dev_lag) %>%
    arrange(dev_lag) %>%
    unique %>%
    .[[1]]

standata_lst <- list(growthmodel_id = 1   # Use weibull rather than loglogistic
                    ,n_data         = usedata_tbl %>% nrow
                    ,n_time         = usedata_tbl %>% select(dev_lag)  %>% unique %>% nrow
                    ,n_cohort       = usedata_tbl %>% select(acc_year) %>% unique %>% nrow
                    ,cohort_id      = get_character_index(usedata_tbl$acc_year)
                    ,cohort_maxtime = cohort_maxtime
                    ,t_value        = t_values
                    ,t_idx          = get_character_index(usedata_tbl$dev_lag)
                    ,premium        = cohort_premium
                    ,loss           = usedata_tbl$cum_loss
                    )

The full Stan file is shown below:

functions {
  real growth_factor_weibull(real t, real omega, real theta) {
    real factor;

    factor = 1 - exp(-(t/theta)^omega);

    return(factor);
  }

  real growth_factor_loglogistic(real t, real omega, real theta) {
    real factor;

    factor = ((t^omega) / (t^omega + theta^omega));

    return(factor);
  }
}

data {
  int<lower=0,upper=1> growthmodel_id;

  int n_data;
  int n_time;
  int n_cohort;

  int cohort_id[n_data];
  int t_idx[n_data];

  int cohort_maxtime[n_cohort];

  real<lower=0> t_value[n_time];

  real premium[n_cohort];
  real loss[n_data];
}

parameters {
  real<lower=0> omega;
  real<lower=0> theta;

  real<lower=0> LR[n_cohort];

  real mu_LR;
  real<lower=0> sd_LR;

  real<lower=0> loss_sd;
}

transformed parameters {
  real gf[n_time];
  real loss_mean[n_cohort, n_time];

  for(i in 1:n_time) {
    if(growthmodel_id == 1) {
      gf[i] = growth_factor_weibull    (t_value[i], omega, theta);
    } else {
      gf[i] = growth_factor_loglogistic(t_value[i], omega, theta);
    }
  }

  for(i in 1:n_cohort) {
    for(j in 1:n_time) {
      loss_mean[i,j] = 0;
    }
  }

  for(i in 1:n_data) {
    loss_mean[cohort_id[i], t_idx[i]] = LR[cohort_id[i]] * premium[cohort_id[i]] * gf[t_idx[i]];
  }
}

model {
  mu_LR ~ normal(0, 0.5);
  sd_LR ~ lognormal(0, 0.5);

  LR ~ lognormal(mu_LR, sd_LR);

  loss_sd ~ lognormal(0, 0.7);

  omega ~ lognormal(0, 1);
  theta ~ lognormal(0, 1);

  for(i in 1:n_data) {
    loss[i] ~ normal(loss_mean[cohort_id[i], t_idx[i]], premium[cohort_id[i]] * loss_sd);
  }
}


generated quantities {
  real mu_LR_exp;
  real<lower=0> loss_sample[n_cohort, n_time];
  real<lower=0> loss_prediction[n_cohort, n_time];
  real<lower=0> step_ratio[n_cohort, n_time];

  real<lower=0> ppc_minLR;
  real<lower=0> ppc_maxLR;

  real<lower=0> ppc_EFC;


  for(i in 1:n_cohort) {
    for(j in 1:n_time) {
      loss_sample[i, j] = LR[i] * premium[i] * gf[t_idx[j]];
      step_ratio [i, j] = 1.0;
    }
  }

  mu_LR_exp = exp(mu_LR);

  for(i in 1:n_data) {
    loss_prediction[cohort_id[i], t_idx[i]] = loss[i];
  }

  for(i in 1:n_cohort) {
    for(j in 2:n_time) {
      step_ratio[i, j] = gf[t_idx[j]] / gf[t_idx[j-1]];
    }
  }

  for(i in 1:n_cohort) {
    for(j in (cohort_maxtime[i]+1):n_time) {
      loss_prediction[i,j] = loss_prediction[i,j-1] * step_ratio[i,j];
    }
  }


  // Create PPC distributions for the max/min of LR
  {
    real tempmin[n_cohort];
    real tempmax[n_cohort];

    for(i in 1:n_cohort) {
      tempmin[i] = LR[i];
      tempmax[i] = LR[i];
    }

    ppc_minLR = min(tempmin);
    ppc_maxLR = max(tempmax);
  }


  // Create total reserve PPC
  ppc_EFC = 0;

  for(i in 1:n_cohort) {
    ppc_EFC = ppc_EFC + loss_prediction[i, n_time] - loss_prediction[i, cohort_maxtime[i]];
  }
}

Need to discuss a number of issues with the above file:

  • Use of functions to select between the Weibull and Loglogistic
  • Use of lognormal to ensure positive values on LR and stddevs
  • Use of normal, not lognormal, for mu_LR as it feeds into a lognormal
  • Use of generated quantities
lc_1_stanfit <- stan(file   = 'losscurves_single.stan'
                    ,data   = standata_lst
                    ,iter   = 500
                    ,chains = 8
                    )

The Stan sample contains no divergent transitions, so that is a good start.

It is always worth checking convergence of the model by checking the \(\hat{R}\) and ensuring it is less than about 1.1

# Plot of convergence statistics
lc_1_draws       <- extract(lc_1_stanfit, permuted = FALSE, inc_warmup = TRUE)
lc_1_monitor_tbl <- as.data.frame(monitor(lc_1_draws, print = FALSE))
lc_1_monitor_tbl <- lc_1_monitor_tbl %>%
    mutate(variable  = rownames(lc_1_monitor_tbl)
          ,parameter = gsub("\\[.*]", "", variable)
           )

ggplot(lc_1_monitor_tbl) +
    aes(x = parameter, y = Rhat, color = parameter) +
    geom_jitter(height = 0, width = 0.2, show.legend = FALSE) +
    geom_hline(aes(yintercept = 1), size = 0.5) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    ylab(expression(hat(italic(R))))
## Warning: Removed 110 rows containing missing values (geom_point).

The \(\hat{R}\) values for the parameter appear to be in or around 1, so that is encouraging.

Finally, we look at some of the n_eff variable also

ggplot(lc_1_monitor_tbl) +
    aes(x = parameter, y = n_eff, color = parameter) +
    geom_jitter(height = 0, width = 0.2, show.legend = FALSE) +
    expand_limits(y = 0) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    xlab("Parameter") +
    ylab(paste0("Effective Sample Count (n_eff)"))

param_root <- c("omega", "theta", "LR", "gf", "loss_sd")

use_vars <- lc_1_monitor_tbl %>%
    filter(parameter %in% param_root) %>%
    .[["variable"]]

Check the traces. There are a large number of parameters in this model fit, so we break them up into groups. First we look at omega, theta and LR

rstan::traceplot(lc_1_stanfit, pars = c("omega", "theta", "LR")) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

Now we look at the traces for gf and loss_sd.

rstan::traceplot(lc_1_stanfit, pars = c("gf", "loss_sd")) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

Now we check the 50% credibility intervals of the parameters

plotdata_tbl <- lc_1_monitor_tbl %>%
    filter(variable %in% use_vars) %>%
    select(mean, `25%`, `50%`, `75%`) %>%
    mutate(variable = factor(use_vars, levels = use_vars))

ggplot(plotdata_tbl) +
    geom_point(aes(x = variable, y = mean)) +
    geom_errorbar(aes(x = variable, ymin = `25%`, ymax = `75%`), width = 0) +
    expand_limits(y = 0) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    xlab("Parameter") +
    ylab("Value")

Now that we have our fit we can start looking at some plots for it. First we look at some very simple sanity-check plots. We look at the full development of an accounting year and see how our fit is doing at tracking that.

fitted_curves_tbl <- extract(lc_1_stanfit)$loss_sample[,1,] %>%
    as_data_frame() %>%
    mutate(iter = 1:n()) %>%
    gather("timelbl", "value", -iter) %>%
    mutate(time = gsub("V", "", timelbl) %>% as.numeric())

ggplot(snapshot_tbl %>% filter(acc_year == 1988)) +
    geom_line (aes(x = time, y = value, group = iter)
              ,data = fitted_curves_tbl, alpha = 0.01) +
    geom_line (aes(x = dev_lag, y = cum_loss), colour = 'red') +
    geom_point(aes(x = dev_lag, y = cum_loss), colour = 'blue') +
    expand_limits(y = 0) +
    scale_y_continuous(labels = scales::dollar) +
    ggtitle("Plot of 1988 Year Loss Development Against Posterior Distribution") +
    xlab("Time") +
    ylab("Loss")

Seems to be doing a reasonable job. Now we use it to predict forward. We start with 1993, where we have development of five years.

predict_cone_tbl <- extract(lc_1_stanfit)$loss_prediction[,6,] %>%
    as_data_frame() %>%
    mutate(iter = 1:n()) %>%
    gather("timelbl", "value", -iter) %>%
    mutate(time = gsub("V", "", timelbl) %>% as.numeric())

plot_predict <- ggplot(carrier_full_tbl %>% filter(grcode == 43, acc_year == '1993')) +
    geom_line (aes(x = time, y = value, group = iter)
              ,data = predict_cone_tbl, alpha = 0.01) +
    geom_line (aes(x = dev_lag, y = cum_loss), colour = 'red') +
    geom_point(aes(x = dev_lag, y = cum_loss), colour = 'blue') +
    expand_limits(y = 0) +
    scale_y_continuous(labels = scales::dollar) +
    ggtitle("Plot of 1993 Year Loss Prediction") +
    xlab("Time") +
    ylab("Loss")

plot(plot_predict)

We now look at a later year, with less claim development.

predict_cone_tbl <- extract(lc_1_stanfit)$loss_prediction[,8,] %>%
    as_data_frame() %>%
    mutate(iter = 1:n()) %>%
    gather("timelbl", "value", -iter) %>%
    mutate(time = gsub("V", "", timelbl) %>% as.numeric())

plot_predict <- ggplot(carrier_full_tbl %>% filter(grcode == 43, acc_year == '1995')) +
    geom_line (aes(x = time, y = value, group = iter)
              ,data = predict_cone_tbl, alpha = 0.01) +
    geom_line (aes(x = dev_lag, y = cum_loss), colour = 'red') +
    geom_point(aes(x = dev_lag, y = cum_loss), colour = 'blue') +
    expand_limits(y = 0) +
    scale_y_continuous(labels = scales::dollar) +
    ggtitle("Plot of 1995 Year Loss Prediction") +
    xlab("Time") +
    ylab("Loss")

plot(plot_predict)

4 Posterior Predictive Checks

A very important part of all this is getting a sense of the aspects of the data that the model is not capturing well. A recommended method for doing this is creating posterior predictive checks (PPCs), that is, assessing the validity of the model in a certain area by comparing the generated values in the sample against the observed values in the dataset.

There are no standard methods for creating PPCs, instead we need to think of different aspects of our data and see how well our model is doing at modelling those idiosyncracies in the data.

We will look at a number of PPCs for the single insurer dataset here.

4.1 Range of Loss Ratios

We first investigate how well the model is doing at capturing the range of Loss Ratios observed in the data: are the largest and smallest predicted Loss Ratios in the model reflecting what we see in the data?

To do this we do a few things: we first add some calculations to the generated quantities block in the Stan file. These calculated values are then compared to what we observed in the data to see how well our model does.

ppc_min_lr <- extract(lc_1_stanfit)$ppc_minLR
ppc_max_lr <- extract(lc_1_stanfit)$ppc_maxLR

lr_tbl <- carrier_full_tbl %>%
    filter(grcode == use_grcode[1]
          ,dev_lag == 10) %>%
    summarise(min_lr = min(loss_ratio)
             ,max_lr = max(loss_ratio))

min_plot <- ggplot() +
    geom_density(aes(x = ppc_min_lr)) +
    geom_vline(aes(xintercept = lr_tbl$min_lr), colour = 'red') +
    xlab("Minimum Loss Ratio") +
    ylab("Probability Density")

max_plot <- ggplot() +
    geom_density(aes(x = ppc_max_lr)) +
    geom_vline(aes(xintercept = lr_tbl$max_lr), colour = 'red') +
    xlab("Maximum Loss Ratio") +
    ylab("Probability Density")


plot_grid(min_plot, max_plot, nrow = 2)

Looking at the above plots, we see that the model is doing reasonably well at capturing the spreads of Loss Ratios.

This is not hugely surprising though, as we have only ten accounting years in the dataset and have seen a decent spread of those already in the dataset.

4.2 Aggregate Reserve Amounts

While it is useful to break down the loss curves into these different cohorts and model each curve separately, from the point of view of an insurance company we care much less about each year’s estimates as we do about the overall amount of money we need to hold back to cover claims.

This presents us with a way to run a check: how well does our model do at estimating the reserves required for all the accounting years put together. We hope that while each accounting year will have mistakes, over-estimates and under-estimates will tend to cancel each other somewhat. Thus, we calculate the total future claims for the book as a whole at 1998 and then compare that to the actual final amounts observed in the data.

As this process may still be a little vague, we will be explicit:

  • For each accounting year \(y\), we look at the current claim level at the point of modelling, 1998. This gives us ten values as we have ten accounting years.
  • We add these ten numbers to calculate \(TCKC_{1998}\), the total of current known claims as at 1998.
  • For each iteration in the sample, we project forward the estimate for the final amount of claims for each accounting year. Summing across the accounting year we end up with a sample of expected final claims, \(EFC_{1998}\).
  • With a sample of this value, we then compare it to the actual, observed values of the variable, \(AFC_{1998}\) and see how the sample values are distributed around the data-calculated value.
tckc <- carrier_snapshot_tbl %>%
    filter(grcode == use_grcode[1]) %>%
    group_by(acc_year) %>%
    filter(dev_lag == max(dev_lag)) %>%
    .[["cum_loss"]] %>%
    sum

afc <- carrier_full_tbl %>%
    filter(grcode == use_grcode[1]) %>%
    group_by(acc_year) %>%
    filter(dev_lag == max(dev_lag)) %>%
    .[["cum_loss"]] %>%
    sum

future_claims <- afc - tckc


ggplot() +
    geom_density(aes(x = extract(lc_1_stanfit)$ppc_EFC)) +
    geom_vline(aes(xintercept = future_claims), colour = 'red') +
    scale_x_continuous(labels = scales::dollar) +
    xlab("Future Claims (,000s)") +
    ylab("Probability Density")

5 Multiple Insurers

stanfile <- 'losscurves_multiple.stan'

ppauto_dt <- ppauto_tbl
setDT(ppauto_dt)

grcodes    <- ppauto_dt[, .(pos_prem = all(earnedpremdir_b > 0)), by = grcode][pos_prem == TRUE, grcode]
use_grcode <- grcodes[1:15]

multi_dt <- ppauto_dt[developmentyear < 1998 &
                      grcode %in% use_grcode, .(grcode    = grcode
                                               ,accyear   = accidentyear
                                               ,premium   = earnedpremdir_b
                                               ,devlag    = developmentlag
                                               ,cumloss   = cumpaidloss_b
                                               ,lossratio = cumpaidloss_b / earnedpremdir_b
                                                )]

cohort_dt <- multi_dt[, .(maxtime = max(devlag), premium = unique(premium)), by = .(grcode, accyear)]
cohort_dt[, cohort_id := .I]
##      grcode accyear maxtime  premium cohort_id
##   1:     43    1988      10      957         1
##   2:     43    1989       9     3695         2
##   3:     43    1990       8     6138         3
##   4:     43    1991       7    17533         4
##   5:     43    1992       6    29341         5
##  ---                                          
## 146:   1767    1993       5 12533746       146
## 147:   1767    1994       4 13590797       147
## 148:   1767    1995       3 14401255       148
## 149:   1767    1996       2 14900682       149
## 150:   1767    1997       1 15065713       150
lst_standata <- list(growthmodel_id = 1   # Use weibull
                    ,n_data         = multi_dt[, .N]
                    ,n_time         = multi_dt[, length(unique(devlag))]
                    ,n_cohort       = cohort_dt[, length(unique(accyear))]
                    ,n_org          = cohort_dt[, length(unique(grcode))]
                    ,n_cohortdata   = cohort_dt[, .N]
                    ,cohort_id      = match(multi_dt$accyear, unique(cohort_dt$accyear))
                    ,org_id         = match(multi_dt$grcode, unique(cohort_dt$grcode))
                    ,t_value        = multi_dt[, sort(unique(devlag))]
                    ,t_idx          = multi_dt[, match(devlag, sort(unique(devlag)))]
                    ,premium        = multi_dt$premium
                    ,loss           = multi_dt$cumloss
                    ,cohort_maxtime = cohort_dt$maxtime
                    )

mi_2_stanmodel <- stan_model(stanfile, verbose = FALSE)

mi_2_stanvb  <- vb(mi_2_stanmodel
                  ,data = lst_standata
                   )
## 
## This is Automatic Differentiation Variational Inference.
## 
## (EXPERIMENTAL ALGORITHM: expect frequent updates to the procedure.)
## 
## Gradient evaluation took 0.000394 seconds
## 1000 iterations under these settings should take 0.394 seconds.
## Adjust your expectations accordingly!
## 
## Begin eta adaptation.
## Iteration:   1 / 250 [  0%]  (Adaptation)
## Iteration:  50 / 250 [ 20%]  (Adaptation)
## Iteration: 100 / 250 [ 40%]  (Adaptation)
## Iteration: 150 / 250 [ 60%]  (Adaptation)
## Iteration: 200 / 250 [ 80%]  (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## 
## Begin stochastic gradient ascent.
##   iter       ELBO   delta_ELBO_mean   delta_ELBO_med   notes 
##    100     -2e+06             1.000            1.000
##    200     -2e+05             3.950            6.901
##    300     -4e+05             2.819            1.000
##    400     -8e+05             2.218            1.000
##    500     -9e+04             3.371            1.000
##    600     -7e+07             2.976            1.000
##    700     -2e+06             8.724            1.000
##    800     -2e+07             7.747            1.000
##    900     -1e+04           196.301            1.000
##   1000     -4e+04           176.744            1.000
##   1100     -9e+03           177.010            3.660   MAY BE DIVERGING... INSPECT ELBO
##   1200     -8e+03           176.327            0.999   MAY BE DIVERGING... INSPECT ELBO
##   1300     -8e+03           176.277            0.999   MAY BE DIVERGING... INSPECT ELBO
##   1400     -8e+03           176.238            0.999   MAY BE DIVERGING... INSPECT ELBO
##   1500     -7e+03           175.442            0.912   MAY BE DIVERGING... INSPECT ELBO
##   1600     -7e+03           175.343            0.729   MAY BE DIVERGING... INSPECT ELBO
##   1700     -7e+03           171.024            0.074   MAY BE DIVERGING... INSPECT ELBO
##   1800     -7e+03           170.933            0.060   MAY BE DIVERGING... INSPECT ELBO
##   1900     -7e+03             0.460            0.025
##   2000     -7e+03             0.388            0.021
##   2100     -7e+03             0.022            0.014
##   2200     -7e+03             0.015            0.013
##   2300     -7e+03             0.009            0.006   MEAN ELBO CONVERGED   MEDIAN ELBO CONVERGED
## 
## Drawing a sample of size 1000 from the approximate posterior... 
## COMPLETED.
mi_2_stanfit <- sampling(mi_2_stanmodel
                        ,data    = lst_standata
                        ,iter    = 500
                        ,chains  = 8
                        ,verbose = FALSE
                         )
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 1).
## 
## Chain 1, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 1, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 1, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 1, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 1, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 1, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 1, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 1, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 1, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 1, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 1, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 1, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 19.9032 seconds (Warm-up)
##                2.20648 seconds (Sampling)
##                22.1097 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 1
##                                                                                          count
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!      2
## Exception thrown at line 101: lognormal_log: Scale parameter is inf, but must be finite!     1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 2).
## 
## Chain 2, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 2, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 2, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 2, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 2, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 2, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 2, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 2, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 2, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 2, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 2, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 2, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 17.4671 seconds (Warm-up)
##                1.96392 seconds (Sampling)
##                19.4311 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 2
##                                                                                          count
## Exception thrown at line 122: lognormal_log: Scale parameter is inf, but must be finite!     3
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!      2
## Exception thrown at line 96: normal_log: Scale parameter is 0, but must be > 0!              1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 3).
## 
## Chain 3, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 3, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 3, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 3, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 3, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 3, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 3, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 3, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 3, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 3, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 3, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 3, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 14.1023 seconds (Warm-up)
##                1.9173 seconds (Sampling)
##                16.0196 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 3
##                                                                                         count
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!     5
## Exception thrown at line 96: normal_log: Scale parameter is 0, but must be > 0!             1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 4).
## 
## Chain 4, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 4, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 4, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 4, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 4, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 4, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 4, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 4, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 4, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 4, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 4, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 4, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 16.251 seconds (Warm-up)
##                2.13724 seconds (Sampling)
##                18.3882 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 4
##                                                                                          count
## Exception thrown at line 122: lognormal_log: Scale parameter is inf, but must be finite!     3
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!      3
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 5).
## 
## Chain 5, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 5, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 5, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 5, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 5, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 5, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 5, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 5, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 5, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 5, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 5, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 5, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 17.4631 seconds (Warm-up)
##                1.83405 seconds (Sampling)
##                19.2971 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 5
##                                                                                          count
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!      3
## Exception thrown at line 122: lognormal_log: Scale parameter is inf, but must be finite!     2
## Exception thrown at line 96: normal_log: Scale parameter is 0, but must be > 0!              1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 6).
## 
## Chain 6, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 6, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 6, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 6, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 6, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 6, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 6, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 6, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 6, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 6, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 6, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 6, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 16.0579 seconds (Warm-up)
##                1.78052 seconds (Sampling)
##                17.8384 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 6
##                                                                                          count
## Exception thrown at line 101: lognormal_log: Scale parameter is 0, but must be > 0!          1
## Exception thrown at line 122: lognormal_log: Scale parameter is inf, but must be finite!     1
## Exception thrown at line 96: normal_log: Scale parameter is 0, but must be > 0!              1
## Exception thrown at line 99: lognormal_log: Scale parameter is inf, but must be finite!      1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 7).
## 
## Chain 7, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 7, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 7, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 7, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 7, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 7, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 7, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 7, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 7, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 7, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 7, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 7, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 16.9446 seconds (Warm-up)
##                1.8352 seconds (Sampling)
##                18.7798 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 7
##                                                                                          count
## Exception thrown at line 122: lognormal_log: Scale parameter is inf, but must be finite!     2
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!      2
## Exception thrown at line 101: lognormal_log: Scale parameter is 0, but must be > 0!          1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
## 
## SAMPLING FOR MODEL 'losscurves_multiple' NOW (CHAIN 8).
## 
## Chain 8, Iteration:   1 / 500 [  0%]  (Warmup)
## Chain 8, Iteration:  50 / 500 [ 10%]  (Warmup)
## Chain 8, Iteration: 100 / 500 [ 20%]  (Warmup)
## Chain 8, Iteration: 150 / 500 [ 30%]  (Warmup)
## Chain 8, Iteration: 200 / 500 [ 40%]  (Warmup)
## Chain 8, Iteration: 250 / 500 [ 50%]  (Warmup)
## Chain 8, Iteration: 251 / 500 [ 50%]  (Sampling)
## Chain 8, Iteration: 300 / 500 [ 60%]  (Sampling)
## Chain 8, Iteration: 350 / 500 [ 70%]  (Sampling)
## Chain 8, Iteration: 400 / 500 [ 80%]  (Sampling)
## Chain 8, Iteration: 450 / 500 [ 90%]  (Sampling)
## Chain 8, Iteration: 500 / 500 [100%]  (Sampling)
##  Elapsed Time: 13.9876 seconds (Warm-up)
##                1.87257 seconds (Sampling)
##                15.8601 seconds (Total)
## The following numerical problems occured the indicated number of times on chain 8
##                                                                                         count
## Exception thrown at line 97: lognormal_log: Scale parameter is inf, but must be finite!     5
## Exception thrown at line 96: normal_log: Scale parameter is 0, but must be > 0!             1
## When a numerical problem occurs, the Hamiltonian proposal gets rejected.
## See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected
## If the number in the 'count' column is small, do not ask about this message on stan-users.
mi_2_draws   <- extract(mi_2_stanfit, permuted = FALSE, inc_warmup = TRUE)
mi_2_monitor <- as.data.frame(monitor(mi_2_draws
                                     ,probs = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975)
                                     ,print = FALSE))
mi_2_monitor$Parameter <- as.factor(gsub("\\[.*]", "", rownames(mi_2_monitor)))


plotdata_dt <- mi_2_monitor[, c('mean', '10%', '50%', '90%')]

setDT(plotdata_dt)
plotdata_dt[, variable := factor(rownames(mi_2_monitor)
                                ,levels = rownames(mi_2_monitor))]
##              mean         10%         50%         90%  variable
##    1:     1.44315     1.39546     1.44255     1.49211  omega[1]
##    2:     1.38639     1.28716     1.38383     1.48703  omega[2]
##    3:     1.25865     1.21228     1.25775     1.30811  omega[3]
##    4:     1.29252     1.23044     1.29240     1.35539  omega[4]
##    5:     1.32334     1.25812     1.32152     1.39288  omega[5]
##   ---                                                          
## 1919:     3.99950     4.00000     4.00000     4.00000 t_098[12]
## 1920:     5.00050     5.00000     5.00000     5.00000 t_098[13]
## 1921:     3.17100     3.00000     3.00000     4.00000 t_098[14]
## 1922:     4.97650     5.00000     5.00000     5.00000 t_098[15]
## 1923: -5460.68370 -5476.92533 -5460.60479 -5445.46766      lp__
ggplot(plotdata_dt[grep("hyper|sd_LR|_exp", variable)]) +
    geom_point(aes(x = variable, y = mean)) +
    geom_errorbar(aes(x = variable, ymin = `10%`, ymax = `90%`), width = 0) +
    expand_limits(y = 0) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    xlab("Parameter") +
    ylab("Value")

ggplot(plotdata_dt[grep("mu_LR_exp", variable)]) +
    geom_point(aes(x = variable, y = mean)) +
    geom_errorbar(aes(x = variable, ymin = `10%`, ymax = `90%`), width = 0) +
    expand_limits(y = 0) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
    xlab("Parameter") +
    ylab("Value")

6 References